home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Collection of Tools & Utilities
/
Collection of Tools and Utilities.iso
/
fortran
/
linpkdrv.zip
/
SUD.FOR
< prev
next >
Wrap
Text File
|
1984-04-30
|
13KB
|
440 lines
C MAIN PROGRAM
INTEGER LUNIT
C ALLOW 5000 UNDERFLOWS.
C CALL TRAPS(0,0,5001,0,0) This is doesn't work on PC's!
C
C OUTPUT UNIT NUMBER
C
LUNIT = 6
C
CALL SQRXX(LUNIT)
STOP
END
SUBROUTINE SQRXX(LUNIT)
INTEGER LUNIT
REAL RX(20,10),R(10,10),XROW(10),Z(10,4),YROW(10),S(10)
REAL SCALE,TINY,HUGE,RHO(2),C(10),SMACH
INTEGER N,P,LDX,LDR,LDZ,CASE
LOGICAL NOTWRT,OFLOW,UFLOW
COMMON SCALE,NOTWRT,OFLOW,UFLOW
LDZ = 10
LDX = 20
LDR = 10
NOTWRT = .TRUE.
OFLOW = .FALSE.
UFLOW = .FALSE.
TINY = SMACH(2)
HUGE = SMACH(3)
SCALE = 1.0E0
DO 60 CASE = 1, 3
GO TO (10, 20, 30), CASE
10 CONTINUE
N = 20
P = 10
GO TO 40
20 CONTINUE
N = 10
P = 4
GO TO 40
30 CONTINUE
N = 10
P = 1
40 CONTINUE
CALL SGETRX(RX,LDX,N,P,S)
WRITE (LUNIT,130) CASE,N,P
IF (NOTWRT) GO TO 50
WRITE (LUNIT,160)
CALL SARRAY(RX,LDX,N,P,P,LUNIT)
50 CONTINUE
CALL SDRUD1(R,LDR,RX,LDX,N,P,XROW,C,S,LUNIT)
CALL SDRUD2(R,LDR,RX,LDX,P,XROW,YROW,Z,LDZ,RHO,C,S,LUNIT)
CALL SDRDD(R,LDR,RX,LDX,P,Z,LDZ,XROW,YROW,RHO,C,S,LUNIT)
60 CONTINUE
CASE = 4
N = 10
P = 4
WRITE (LUNIT,130) CASE,N,P
WRITE (LUNIT,140)
CALL SGETRX(RX,LDX,N,P,S)
DO 80 J = 1, P
DO 70 I = 1, N
RX(I,J) = HUGE*RX(I,J)
70 CONTINUE
80 CONTINUE
SCALE = HUGE
OFLOW = .TRUE.
IF (NOTWRT) GO TO 90
WRITE (LUNIT,160)
CALL SARRAY(RX,LDX,N,P,P,LUNIT)
90 CONTINUE
CALL SDRUD1(R,LDR,RX,LDX,N,P,XROW,C,S,LUNIT)
CALL SDRUD2(R,LDR,RX,LDX,P,XROW,YROW,Z,LDZ,RHO,C,S,LUNIT)
CALL SDRDD(R,LDR,RX,LDX,P,Z,LDZ,XROW,YROW,RHO,C,S,LUNIT)
OFLOW = .FALSE.
CASE = 5
N = 10
P = 4
WRITE (LUNIT,130) CASE,N,P
WRITE (LUNIT,150)
CALL SGETRX(RX,LDX,N,P,S)
DO 110 J = 1, P
DO 100 I = 1, N
RX(I,J) = TINY*RX(I,J)
100 CONTINUE
110 CONTINUE
SCALE = TINY
UFLOW = .TRUE.
IF (NOTWRT) GO TO 120
WRITE (LUNIT,160)
CALL SARRAY(RX,LDX,N,P,P,LUNIT)
120 CONTINUE
CALL SDRUD1(R,LDR,RX,LDX,N,P,XROW,C,S,LUNIT)
CALL SDRUD2(R,LDR,RX,LDX,P,XROW,YROW,Z,LDZ,RHO,C,S,LUNIT)
CALL SDRDD(R,LDR,RX,LDX,P,Z,LDZ,XROW,YROW,RHO,C,S,LUNIT)
RETURN
C
130 FORMAT (11H1 CASE =, I2, 5X, 3HN =, I2, 5X, 3HP =, I2 /////)
140 FORMAT (22H OVERFLOW TEST /////)
150 FORMAT (24H UNDERFLOW TEST /////)
160 FORMAT (5H RX)
END
SUBROUTINE SARRAY(A,LDA,M,N,NNL,LUNIT)
INTEGER LDA,M,N,NNL,LUNIT
REAL A(LDA,1)
C X
C FORTRAN IABS,MIN0
C
INTEGER I,J,K,KU,NL
NL = IABS(NNL)
IF (NNL .LT. 0) GO TO 30
DO 20 I = 1, M
WRITE (LUNIT,70)
DO 10 K = 1, N, NL
KU = MIN0(K+NL-1,N)
WRITE (LUNIT,70) (A(I,J), J = K, KU)
10 CONTINUE
20 CONTINUE
GO TO 60
30 CONTINUE
DO 50 J = 1, N
WRITE (LUNIT,70)
DO 40 K = 1, M, NL
KU = MIN0(K+NL-1,M)
WRITE (LUNIT,70) (A(I,J), I = K, KU)
40 CONTINUE
50 CONTINUE
60 CONTINUE
RETURN
C
70 FORMAT (1H , 4(2E13.6, 4X))
END
SUBROUTINE SDRDD(R,LDR,RX,LDX,P,Z,LDZ,XROW,YROW,RHO,C,S,LUNIT)
INTEGER LDR,LDX,LDZ,P,LUNIT
REAL R(LDR,1),RX(LDX,1),Z(LDZ,1),XROW(1),YROW(1),S(1)
REAL RHO(1),SCALE,C(1)
LOGICAL NOTWRT,OFLOW,UFLOW
COMMON SCALE,NOTWRT,OFLOW,UFLOW
INTEGER I,J,JP1,PM1
WRITE (LUNIT,50)
CALL SCHDD(R,LDR,P,XROW,Z,LDZ,2,YROW,RHO,C,S,INFO)
PM1 = P - 1
IF (PM1 .LT. 1) GO TO 30
DO 20 J = 1, PM1
JP1 = J + 1
DO 10 I = JP1, P
R(I,J) = 0.0E0
10 CONTINUE
20 CONTINUE
30 CONTINUE
IF (NOTWRT) GO TO 40
WRITE (LUNIT,80)
CALL SARRAY(R,LDR,P,P,P,LUNIT)
WRITE (LUNIT,60)
CALL SARRAY(Z,LDZ,P,2,2,LUNIT)
WRITE (LUNIT,70) (RHO(I), I = 1, 2)
40 CONTINUE
CALL SMPDD(R,LDR,P,RX,LDX,Z,Z(1,3),LDZ,RHO,LUNIT)
RETURN
50 FORMAT ( /////
* 46H STEP THREE DOWNDATING XROW,YROW AND Z, ///)
60 FORMAT ( /// 4H Z)
70 FORMAT ( /// 6H RHO // 1H , 2E13.6)
80 FORMAT (4H R)
END
SUBROUTINE SDRUD1(R,LDR,RX,LDX,N,P,XROW,C,S,LUNIT)
INTEGER N,P,LDR,LDX,LUNIT
REAL R(LDR,1),RX(LDX,1),XROW(1),S(1)
REAL C(1)
LOGICAL NOTWRT,OFLOW
REAL RELM,XELM,Y,Z
REAL XMAX,RMAX,TEST,T,SCALE,SMACH
COMMON SCALE,NOTWRT,OFLOW,UFLOW
DO 20 I = 1, P
DO 10 J = 1, P
R(I,J) = 0.0E0
10 CONTINUE
20 CONTINUE
DO 40 I = 1, N
DO 30 J = 1, P
XROW(J) = RX(I,J)
30 CONTINUE
CALL SCHUD(R,LDR,P,XROW,Z,1,0,Y,Y,C,S)
40 CONTINUE
WRITE (LUNIT,100)
IF (NOTWRT) GO TO 50
WRITE (LUNIT,120)
CALL SARRAY(R,LDR,P,P,P,LUNIT)
50 CONTINUE
RMAX = 0.0E0
XMAX = 0.0E0
DO 90 I = 1, P
DO 80 J = 1, I
RELM = 0.0E0
XELM = 0.0E0
NIM = MIN0(I,J)
DO 60 K = 1, NIM
RELM = RELM + R(K,I)/SCALE*(R(K,J)/SCALE)
60 CONTINUE
DO 70 K = 1, N
XELM = XELM + RX(K,I)/SCALE*(RX(K,J)/SCALE)
70 CONTINUE
T = AMAX1(ABS(XELM),0.0E0)
XMAX = AMAX1(XMAX,T)
T = AMAX1(ABS(XELM-RELM),0.0E0)
RMAX = AMAX1(RMAX,T)
80 CONTINUE
90 CONTINUE
TEST = RMAX/XMAX/SMACH(1)
WRITE (LUNIT,110) TEST
RETURN
C
100 FORMAT ( ///// 25H STEP ONE UPDATING X ///)
110 FORMAT ( /// 15H STATISTICS //
* 42H RH*R ............................, E10.2)
120 FORMAT (4H R)
END
SUBROUTINE SDRUD2(R,LDR,RX,LDX,P,XROW,YROW,Z,LDZ,RHO,C,S,LUNIT)
INTEGER LDR,LDX,LDZ,P,LUNIT
REAL R(LDR,1),RX(LDX,1),XROW(1),YROW(1),Z(LDZ,1),S(1)
REAL RHO(1),C(1)
REAL SCALE
LOGICAL NOTWRT,OFLOW,UFLOW
COMMON SCALE,NOTWRT,OFLOW,UFLOW
DO 20 I = 1, P
DO 10 J = 1, P
RX(I,J) = R(I,J)
10 CONTINUE
20 CONTINUE
DO 40 I = 1, P
Z(I,1) = 0.0E0
DO 30 J = I, P
Z(I,1) = Z(I,1) + R(I,J)
30 CONTINUE
Z(I,2) = Z(I,1) + SCALE
Z(I,3) = Z(I,1)
Z(I,4) = Z(I,2)
40 CONTINUE
DO 60 I = 1, P
XROW(I) = 0.0E0
DO 50 J = 1, I
XROW(I) = XROW(I) + R(J,I)
50 CONTINUE
60 CONTINUE
YROW(1) = 0.0E0
DO 70 I = 1, P
YROW(1) = YROW(1) + XROW(I)
70 CONTINUE
YROW(2) = YROW(1) - SCALE
RHO(1) = SCALE
RHO(2) = SCALE
IF (NOTWRT) GO TO 80
WRITE (LUNIT,100)
CALL SARRAY(XROW,P,P,1,-P,LUNIT)
WRITE (LUNIT,110)
CALL SARRAY(Z,LDZ,P,2,2,LUNIT)
WRITE (LUNIT,120)
CALL SARRAY(YROW,2,2,1,-2,LUNIT)
80 CONTINUE
CALL SCHUD(R,LDR,P,XROW,Z,LDZ,2,YROW,RHO,C,S)
WRITE (LUNIT,130)
IF (NOTWRT) GO TO 90
WRITE (LUNIT,150)
CALL SARRAY(R,LDR,P,P,P,LUNIT)
WRITE (LUNIT,110)
CALL SARRAY(Z,LDZ,P,2,2,LUNIT)
WRITE (LUNIT,140) (RHO(I), I = 1, 2)
90 CONTINUE
CALL SMPUD(R,LDR,RX,LDX,P,XROW,Z,LDZ,RHO,LUNIT)
RETURN
C
100 FORMAT ( ///// 7H XROW)
110 FORMAT ( /// 4H Z)
120 FORMAT ( /// 7H YROW)
130 FORMAT ( ///// 41H STEP TWO UPDATING XROW, YROW AND Z ///)
140 FORMAT ( /// 6H RHO // 1H , 2E13.6)
150 FORMAT (4H R)
END
SUBROUTINE SGETRX(X,LDX,N,P,S)
INTEGER N,P,LDX
REAL X(LDX,1),S(1)
INTEGER PD2,PD2D1,I
PD2 = MAX0(P/2,1)
PD2D1 =